Gyrokinetic simulation of entropy cascade 
in two-dimensional electrostatic turbulence 
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Two-dimensional electrostatic turbulence in magnetized weakly-collisional plasmas exhibits a cascade of 
entropy in phase space [Phys. Rev. Lett. 103, 015003 (2009)]. At scales smaller than the gyroradius, this cascade 
is characterized by the dimensionless ratio D of the collision time to the eddy turnover time measured at the scale 
of the thermal Larmor radius. When D » 1, a broad spectrum of fluctuations at sub-Larmor scales is found in 
both position and velocity space. The distribution function develops structure as a function of Vj_, the velocity 
coordinate perpendicular to the local magnetic field. The cascade shows a local-scale nonlinear interaction in 
both position and velocity spaces, and Kolmogorov's scaling theory can be extended into phase space. 
Keywords: electrostatic turbulence, gyrokinetic theory, nonlinear phase mixing, kinetic dissipation, decaying 

turbulence simulation, Kelvin-Helmholtz instability 



1. Introduction 

Fluid turbulence can be described mostly by Navier- 
Stokes equations, which consist of partial differential equa- 
tions in the position space — thus in 2D or 3D space. Such 
a description is possible because we may regard the fluid 
to be sufficiently collisional and to be locally in thermo- 
dynamic equilibrium. From the outer driving scales to the 
small dissipative scales, the energy (or the enstrophy in the 
2D case) cascades in the wave-number space |T). The iner- 
tial range dynamics are basically governed by the nonlin- 
ear term. The fluxes are finally dissipated by the viscosity 
described by a diffusion operator, in the dissipation range 
of wave-number space. There is a dimensionless number 
called the Reynolds number, which characterizes the scale 
separation of the turbulent dynamics. The ratio of the dis- 
sipation scale to the outer scale is in general determined by 
a fractional power of the Reynolds number. 

Kinetic turbulence is less well understood. We have 
some knowledge especially from the extensive study of fu- 
sion plasmas |[2}0] and recently of space plasmas J5]^l, 
but not to the same extent as fluid turbulence is known, par- 
ticularly in the sense described above. Since most plasmas 
are in a collisionless or weakly-collisional environment, 
they are not in local thermodynamical equilibrium. In this 
case we have to invoke kinetic description in the phase 
space, which enables us to use Vlasov or Boltzmann the- 
ory, or their magnetized low-frequency limit, the gyroki- 
netics. Dissipation in the kinetic system occurs because 
of collisions. Without collisions, the entropy is conserved 
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and the system is in principle reversible in a thermody- 
namic sense IflOl . Wave resonances or Landau damping do 
not by themselves constitute dissipation since no entropy 
is generated unless we include the effect of collisions. Ir- 
reversibility becomes the key issue here. 

In this context, the following question arises. Col- 
lisions are usually described by a diffusion operator in 
the velocity space IfTTI . How do collisions in a weakly 
collisional system determine the wave-number as well as 
velocity-space cutoff of the inertial range? In this pa- 
per, following Ref. fl~2l . we show detailed numerical evi- 
dence of the coupling between position- and velocity-space 
structures, and the consequent achievement of the colli- 
sional dissipation described by the velocity-space diffu- 
sion operator. We introduce a new dimensionless number 
D, which, analogously to the Reynolds number in Navier- 
Stokes equations, characterizes the scale separation be- 
tween the outer scale and the dissipation scale in the ki- 
netic theory. We emphasize that we are concerned only 
with scales smaller than the Larmor radius. 

In order to focus on the nonlinear interaction, we omit 
Landau damping by ignoring variation along the field line. 
The system retains two collisionless invariants as the 2D 
Navier-Stokes equations do. These two invariants cannot 
share the same local-interaction space in a Kolmogorov- 
like phenomenology, and thus lead to a dual cascade (direct 
and inverse cascades) OTL3TfT6ll . In this paper we focus on 
the direct cascade only and the discussion of the inverse 
cascade is left for future publication. 

This paper is organized as follows. As a reduced ki- 



netic model for plasma turbulence, we first introduce the 
gyrokinetic (GK) equation briefly and describe its basic na- 
ture in Sec. [2] Scaling relations of the turbulent cascade are 
also briefly summarized. Section [3] shows the evidence of 
kinetic turbulent cascade by means of the numerical simu- 
lation of the decaying turbulence. In addition to the repro- 
duction of the theoretical prediction, we show how much 
resolution is needed in both position and velocity space to 
achieve the proper dissipation and show the trend that the 
amount of dissipation is asymptotically independent of the 
collision frequency. We also present the characteristics of 
the nonlinear triad interaction in wave-number space and 
show that the assumption of local-scale interaction is sup- 
ported. We finally summarize our results in Sec. [4] 

2. Gyrokinetics 
2.1 Equations 

We first introduce the GK model briefly Ifl7ti20ll . 
Since we are going to look at the turbulence in the mag- 
netized plasmas, the dynamics of interest is much slower 
than particles' gyromotion. We may average over the gy- 
romotion and the gyroangle can be ignored due to this gy- 
roaveraging. A GK system has 3 spatial coordinates (x, y, 
z), and 2 velocity coordinates (vx, \'\\), where || and ± stand 
for parallel and perpendicular directions to the background 
magnetic field. It is convenient to distinguish between the 
particle coordinate r and the gyrocenter coordinate R IfTTll . 
These coordinates are connected by 



R = r + 



vxz 
CI ' 



(1) 



where z is the unit vector along the background magnetic 
field and Q is the gyrofrequency. 

We further reduce the GK equation into 2D in position 
space, or 4D in phase space, by ignoring variation along 
the mean field (k\\ = 0). This is useful because we not only 
reduce the dimension but also remove the Landau damping 
from the system. The resultant equation is 



+ —{(<P)R,h) = — 7" + (C[h]) R , 

at B r () at 



(2) 



where h is the non-Boltzmann part of the perturbed ion 
distribution function, ip is the electrostatic potential, Bq is 
the background magnetic field aligned with the z-axis, {-)r 
is the gyroaverage holding the guiding center position R 
constant, {/, g] — {z x V/) • Vg is the standard Poisson 
bracket, q is the charge, and Tq is the temperature of the 
background Maxwellian Fq. The collision operator C[h] 
we use in our simulations describes pitch-angle scattering 
and energy diffusion with proper conservation properties 
(see El). 

To calculate the self-consistent potential at these small 
scales, we use the quasineutrality condition 



where (-} r denotes the gyroaverage at fixed particle posi- 
tion r, Q — Yisl^nos/Tos for Boltzmann-response (3D) 
electrons or Q — qfnoi/Toi for no-response (2D) electrons, 
«o is the density corresponding to Fq, and s and i are the 
species indices. In this paper we use no response electrons 
since electrons cannot contribute to the potential if k\\ = 
exactly Il22l . We have also made a separate two-species 
simulation with real mass ratio and obtained almost iden- 
tical results in the regime k ± p e < 1 with negligible density 
fluctuations, where p e is the electron Larmor radius. 
There are two collisionless invariants, 



e = f 2 (1 ~ ro)|( ^ 2 ' 



<p 2 dr, 



(4) 
(5) 



where the subscript k is for the Fourier coefficients, To = 
Io(b)e~ h , Iq is the modified Bessel function, b = k]_p 2 /2, 
and p is the ion thermal Larmor radius. The invariant 3B 
is proportional to the negative of the perturbed part of the 
entropy of the system, - J f In / dr dv ||2][23), where / is 
the full distribution function. Here we will refer to 38 as 
"entropy" to emphasize this connection, although in ther- 
modynamic terms, it is better interpreted as a generalized 
free energy Il24ll25l . The second invariant € implies the 
existence of a dual cascade in the 2D system, which we 
do not discuss in this paper. Note that these quantities are 
conserved only in the collisionless limit. Any decrease of 
28 due to collisions corresponds to the creation of entropy 
and heating, hence the irreversibility. 

2.2 Turbulent scaling 

A scaling theory of the entropy cascade in the sub- 
Larmor scale range it «: p) can be developed in a way 
reminiscent of the Kolmogorov-style turbulence theories. 
From the conservation of SB in the collisionless limit, we 
may assume that the entropy is going to cascade scale by 
scale without dissipation in the inertial range. Note that the 
quantity which cascades forward (to smaller scales) is 38, 
not €. This has been demonstrated by applying Fj0rtoft's 
argument on the dual cascade 11311251 . 

The nonlinear term [the Poisson bracket in <j2j] that 
causes the entropy transfer introduces velocity-space struc- 
ture simultaneously. For small-scale electric fields, parti- 
cles with different gyroradii execute different E x B mo- 
tions because they see different effective potentials; this 
leads to nonlinear phase mixing and other novel phenom- 
ena 12211261 . It can be argued that the nonlinear phase mix- 
ing produces the following relationship between position- 
space and velocity-space scales: 



6v ± e 

Vth P ' 



(6) 



Q<p = q \ (h) r dv 



(3) 



Dimensional arguments then lead to the following 
scalings of the spectra (see Refs. |[T2]|25]|27| for details), 



Table 1 Estimated dimensionless number D in various physical 
systems. Ion scale turbulence at LAPD |28), DIII-D (29) 
and TFTR (30), ETG turbulence at NSTX QT), and solar 
wind observation (5). Eddy turnover time t for the ETG 
case is evaluated from simulation data 1321 . 





v [sec ] 


t [sec] 


D 


LAPD 


2.1 ■ 10 6 


> 5.0 


• 10~ 6 


<0.1 


DIII-D 


3300 


7.5 ■ 


10~ 6 


40 


TFTR 


120 


1.5 ■ 


io- 5 


500 


ETG @ NSTX 


10,800 


3.5 ■ 


io- 7 


270 


Solar wind 


2.3 ■ 10~ 7 


190 


2.3 ■ 10 4 



E h (k ± )Kk- 4 ' 3 , Ev(k ± )Kkl LU ", 



10/3 



where 



\k ± \=k± 



(7) 

(8) 
(9) 



Note that the total entropy © can be expressed as 36 = 
/ [E h (kJ - E v (kj_)] dk ± . 

At small position-space scales, collisions tend to dom- 
inate as velocity-space structure becomes finer due to ©. 
It can be shown [12, 25 . 27 1 that the cutoffs are 



6v ±c 1 



Vth k ±c p 



D 



-3/5 



where 



D 



1 

VT n 



(10) 



(11) 



the quantity v is the ion collision frequency and t p is a 
turnover time at the scale p. 

We introduced a dimensionless number D which rep- 
resents the ratio between the outer scale, in this case the 
thermal Larmor radius, and the dissipation scale valid for 
the kinetic theory. This is simply the ratio of the colli- 
sion time and the eddy turnover time at the ion gyroscale. 
It describes how many times eddies at the Larmor radius 
scale turn over before the cumulative effect of small-angle 
collisions smears out the resulting structures in velocity 
space at this scale. Analogous to the Reynolds number, 
D becomes large when dissipation is weak, i.e., when the 
collision frequency is small. We emphasize here that this 
number is amplitude dependent since turnover time is de- 
termined by the E x B velocity where £ is a fluctuation. 

Table \T\ shows D values in more realistic and geomet- 
rically complicated physical systems. The first three are 
the ion-scale turbulence in the basic plasma device LAPD 
It28l and in the fusion devices DIII-D J29) and TFTR QUI 



The fourth one is the electron-temperature-gradient (ETG) 
driven turbulence in NSTX, where v is estimated from ex- 
perimental conditions PD while t is from the streamer- 
dominated results of simulation 11321 at k ± p e = 0.2. The 
last one is for the solar wind turbulence for which the data 
from Cluster spacecraft at 19 earth radii is used Q. In 
the solar wind, the collision frequency is extremely small 
because of the small density, leading to extremely large D. 

3. Numerical simulation 
3.1 Code: AstroGK 

AstroGK is an Eulerian initial value solver for the GK 
equations in 5D phase space in general, although we only 
use 4 dimensions in this study. Namely, we discretize the 
velocity space into grids as standard fluid codes do in the 
position space. It evolves the function 



= h- 



To 



(12) 



according to (f2]i and ©. In order to make a 4D simula- 
tion using the 5D code, we used only 3 grid points along 
the field line (in z) with a very long box size and a peri- 
odic boundary condition. Starting from the homogeneous 
initial condition along z, we confirmed that no structure 
developed in this direction. 

The code uses a Fourier pseudo-spectral scheme for 
real space dimensions perpendicular to the mean mag- 
netic field and Legendre collocation scheme in the two- 
dimensional velocity space integration. The velocity-space 
grid is taken with respect to particle energy s = v 2 and 
A = v\le, which makes our grid points radially distributed 
in the vy-Vj. plane 11331 . For the collision operator we use a 
conservative finite difference scheme 041 . 

Time integration is made using the 3rd order Adams- 
Bashforth scheme for the nonlinear term. The linearized 
collision operator is treated by the first-order implicit Euler 
scheme with Sherman-Morrison formula for the moments- 
conserving corrections I34II35I . 

AstroGK is written in Fortran 95 and is parallelized 
using MPI. Our biggest run (D2 in Table [2]i cost about 36 
wall-clock hours using 8,192 processor cores. 

3.2 Initial condition 

Here we describe the setup of the direct cascade sim- 
ulations. We use the simple straight slab geometry in a 
homogeneous background in the box size L x = L y = 27rp. 
A series of decaying turbulence simulations were carried 
out, with the initial condition put in \k x p\, \k v p\ = 2 scale: 



ginit - g0 



2x 2y 

COS 1- COS >r X(X,y) 

P P 



F . 



(13) 



where Fo is a Maxwellian. Here, go is a constant and 
X(x,y) is small-amplitude white noise superimposed on all 
Fourier harmonics. From (O and d 121) . we can calculate 
i^init- The initial condition ( fT"3l is unstable to the kinetic 



Table 2 Index of the runs described in Sec. [3] Dimensionless 
number D is given for the initial stage of developed turbu- 
lence (t/T mit = 10), and cutoff wave number is estimated 
by (16) (see Sec.[33]and Fig.[3). 



Run 


VT init 




D 


k ±c p 


N X XNy 


N e x 2N A 


A 


9.3 • 10" 


-3 


32 


16 


64 2 


32 2 


B 


5.6-10 


-3 


48 


20 


64 2 


32 2 


CI 


1.9-10 


-3 


118 


35 


128 2 


16 2 


C2 


1.9-10 


-3 


118 


35 


128 2 


64 2 


C3 


1.9-10 


-3 


118 


35 


128 2 


128 2 


Dl 


7.4- 10 


-4 


440 


77 


128 2 


128 2 


D2 


7.4 • 10" 


-4 


440 


77 


256 2 


128 2 



version of the Kelvin-Helmholtz (KH) instability, which 
self-consistently evolves into turbulence. It is a completely 
autonomous system and there is no entropy input during 
the time evolution. We emphasize that there is no non- 
Maxwellian velocity-space structure in the initial condi- 
tion. All fine structure later observed in the velocity space 
is created from nonlinear interaction described in Sec. [2] 

The results reported below were obtained in the series 
of runs indexed in Table [2] where x N y is number of 
collocation points in the position space and N e x 2N A is 
the number of grid points in velocity space — the factor 
of 2 corresponds to the sign of vy = ± y/s(l - A). For each 
collision frequency, D is calculated using the turnover time 
taken from the simulation (see Sec. 13.31 and Fig. [3]). The 
cutoff scale is listed for each case here, and we determine 
the grid resolution so that it resolves the smallest position- 
and velocity-space structures. 

3.3 Time evolution 

Time evolution of the collisionless invariants and the 
amplitude of several lower Fourier modes |(^| are shown 
in Figs. QJO where T; n i t is the initial turnover time of the 
gyroaveraged E x B motion defined by 

27rB 



Tinit — 



ck x k y \\{if in 
with \k x p\, \k y p\ = 2 and 



>*ll 



\\(<P)r\\ = 



"o II 



\(<p) R \ 2 F dvdR 



1/2 



(14) 



(15) 



In Fig.Q] the decrease of 38 corresponds to the creation of 
entropy since 38 is proportional to the negative of the per- 
turbed entropy. Initially, while instability grows (t /Tinit ^ 
8, see also Fig. 13, 38 decays slowly with a rate proportional 
to the collision frequency v. Then as the instability satu- 
rates around f /Ti n it - 9, 38 starts to decay very rapidly. This 
rapid decay of 38 suggests nonlinear creation of small scale 
structures in velocity space as well as in position space. 
Some time after that (f/Ti n it 20), 38 starts to decay expo- 
nentially suggesting the end of turbulence. The constancy 
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Fig. 1 



Time evolution of (a) 38 and € and J5J], and (b) 
rate of change of 38, normalized to initial 38. Labels A- 
D2 correspond to runs indexed in Table [2] Evolution of 
£ does not differ among runs significantly, and is repre- 
sented by the run D2. 



10' 



« io-; 

^ 10 "; 

9: 10! 

10 
10/ 
10 



10 u 

10". 

io-; 

10- 



10 































































































k 


,=o 




1 












h 




sT 


o 








1=0 


y~ 


2 - 








k y = 
k - 




























2 4 6 8 10 12 14 16 18 20 



t/-. 



init 



Fig. 2 Time evolution of several lower Fourier modes ip^ from 
the run D2 (See Table|2j. 



of € implies the presence of an inverse cascade 1271 . which 
will be discussed elsewhere. 

Time evolution of the rate of change of 38 is shown 
in Fig. Hlb)- Note that the decay of 38 is accounted for 
perfectly by collisional entropy production — i.e. entropy 
balance J51 is satisfied. This is guaranteed under the con- 
trolled velocity-space numerical scheme ll33l . In all runs 
except for the one with the smallest D (run A in Table |2j 
D = 32), the rate of change of 38 is almost identical in the 
nonlinear stage (f/Tinit <: 9). The slight horizontal shift of 
the lines is due to the different timing of linear mode satu- 
ration because of the random-noise admixture in the initial 
conditions. This plot suggests that the asymptotic amount 
of dissipation is finite as collision frequency tends to zero, 
so the dissipation rate tends to a value asymptotically in- 
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Fig. 3 Time evolution of the dimensionless number D [see 
dl 1H. Labels B-D2 correspond to runs indexed in Table 
[2] Values of D in Table[2]is evaluated at t/r init = 10. 



dependent of the collisionality, but rather determined by 
the nonlinear dynamics of the turbulence. Also, the falling 
phase of dlB/dt corresponds to the development of the tur- 
bulence, and we may say that the turbulence has fully de- 
veloped by f/Tinit = 10. In the following sections, we will 
concentrate on the asymptotic cases, which have similar 
rates of change of 38 (runs B-D). 

The dimensionless number D is defined as follows. 
As is seen in Fig. [2j the simple use of a single represen- 
tative wave number k ± p = 2 yields significant fluctuation 
of D which is not necessarily physical. On the other hand, 
the k ± p = 1 component is created by the KH instability 
as a part of the inverse cascade, and is also irrelevant for 
characterizing direct cascade. Thus we introduced a def- 
inition of the turnover time at scale I - np, denoted by 
T np , which uses all k ± p > 2 components of the fluctuating 
potential. Using r ffp , the time evolution of D is shown in 
Fig. [3] The values are shown for three resolved cases from 
f/Ti n it ^ 8 since the turbulent spectra start to develop after 
that. At f/Ti n it = 10 the dimensionless number D takes val- 
ues shown for each run in Table[2] The cutoff wave number 
is then obtained from 
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Fig. 4 Time evolution of the wave-number spectra of 
from the run D2 (see TableO. 
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Fig. 5 Time-averaged normalized wave-number (Fourier) spec- 
tra (a) E h (k x )/W and (b) ^(yfcJ/38 [cf. Q] for the runs 
indexed in Table [2] Theoretically predicted slopes are 
given for comparison. 



k ±c p = aD 3 ' 5 . 



(16) 



The value a — 2 corresponds to the particular initial con- 
dition and domain size used here. 

3.4 Wave-number spectra 

We first show the time evolution of the wave-number 
spectra of E/,(k ± ) averaged over wave-number shell [de- 
fined by ©] in Fig. [4] At f/T; n it = 9 the spectrum shows an 
overall steeper slope as it has not yet entered the fully de- 
veloped turbulence stage. From f/Tinit = 1 1 to 14 the spec- 
tra show a self-similar shape consistent with the inertial- 
range scaling (0. The spectra at f/T; n i t = 17 and 20 clearly 
have an inertial range; the steeper falloff at large k ± sug- 
gests a dissipation range. Note that the implied cutoff scale 
is consistent with the estimate from the time-evolving di- 
mensionless number D [see Fig.[3]and (U6H , 

Figure \5\ shows the wave-number spectra of Eh{ki_) 



and E^ik^) [see ([8]) and (0] of the asymptotic runs (B- 
D in Table |2]i which are averaged over wave-number shell 
and normalized by the value of total 38 at each time, and 
then averaged over time 10 < f/T; n it < 15. It is hard to 
identify the dissipation range as the time average is taken 
over the period when the cutoff is near the highest k ± end 
so that inertial range can be taken as wide as possible. Un- 
clear dissipation range in the time-averaged spectra is ac- 
ceptable as exponential fall-off in this regime implies ex- 
ponentially small dissipation apart from the cutoff. Three 
resolved cases B, C2 and D2 lie on top of one another and 
are consistent with the theoretical slopes (£/, oc A:~ 4 ^ 3 and 
oc A;^ 10 ^ 3 ). In contrast, the under-resolved runs in the 
velocity space (CI) and in the position space (Dl) have 
shallower slopes than others. With the increase of D the 
dissipation cutoff extends to smaller scales, thus higher res- 




50 100 



Fig. 6 Time-averaged normalized 2D spectrum E g (k ± ,p) [de- 
fined by J17H taken from the run D2 (See Table[2j- 



olution is required both in position and velocity space. The 
simulation shows a consistent result with the scaling theory 
including the dissipation cutoff. 

3.5 Velocity-space spectra 

In order to characterize the velocity-space structure 
quantitatively, Plunk et al. Il27l proposed a dual space vari- 
able p to the velocity variable v ± and introduced a 2D spec- 
trum in (k ± ,p) space: 



E g (k ± ,p)= J] P\UP)\ 2 

\k ± \=k x 



(17) 



where gk(p) = J Jo(pv±)gk(v) dv is the Hankel transform. 

Figure [6] shows the time average of the normalized 
2D spectrum E g (k ± ,p) taken from the run D2 (see Table 
E). We first note that there is a gap in the low-k ± , high- 
p region. This gap appears because high-/? component is 
not easily created with a low-k ± mode, while the high-A: ± 
component can be created even for low-/? as Navier-Stokes 
turbulence does with Maxwellian velocity profile. This is 
likely due to the fact that nonlinear phase mixing acts more 
strongly at sub-Larmor scales. The spectral contour devel- 
ops along the diagonal pv t ^ = k ± p, which reflects the fact 
that the correlation of position- and velocity-space struc- 
ture follows our conjecture ((6) very well. Thus, compa- 
rable resolution in velocity and position space is required, 
especially when we investigate smaller scales than the Lar- 
mor radius. 

In Ref. l27l the asymptotic spectra are derived for 
k_i_p » pv t h and k ± p <& pv t h limits, namely 



fcV 1/3 (k^p» P v th ) 
E g (k ± ,p) cc\ 

5 ^k ± (k ± p « pvth) 



(18) 



Figure|7]shows various (a) horizontal and (b) vertical slices 
of Fig. [6] Depending on the value of p in Fig. [7ja), we ob- 
serve a slope approaching k^ in the high-£ ± regime. On 
the other hand, the £~^ 3 slope for k ± p <K pv± seems diffi- 
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Fig. 7 Averaged asymptotic spectra from the run D2 (See Table 
[2J ( a ) as a function of k ± for various fixed p and (b) as a 
function of p for various fixed k±. 
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Fig. 8 Time-averaged normalized velocity-space (Hankel) 
spectrum E g (p) = J E g (k ± ,p)dk L [see also Q7)] for the 
runs in Table [2] 



cult to obtain, reflecting the gap in the low-£ ± , high-/? re- 
gion seen in Fig. [6] Figure|7jb) shows a rather good slope 
consistent with that predicted for the low-/? limit ( oc /?~'/ 3 ) 
even for the smallest k ± in the figure. The high-/? slope 
is steeper than the prediction (fT8l for the k ± p = 10 case 
because of the remnant of the gap in Fig. [6] The agree- 
ment with the asymptotic expectation becomes good for 
k±p k. 30 even though we don't have a wide asymptotic 
range in this region. The overall spectra show a reason- 
able symmetry in k± and p except for the structure coming 
from the gap, which indicates consistency with the theoret- 
ical prediction (fT8l . 

Figure [8] shows the time-averaged spectra in the Han- 
kel space E g {p) = J E g (k ± ,p) dk ± normalized to 3B at each 
time. The theoretical expectation is Il27l 



E g (p) ' 



,-4/3 



(19) 



which logically follows from the first spectrum in (01 
and pv t h ~ k x p [equivalent to ©]. The numerical re- 
sult again shows approximate consistency with the theo- 
retical prediction and confirms that small-scale structure 
is formed in the velocity space. The large hump in the 
low-/? regime (p < 3) is due to long-wave-length modes, 
which have significantly larger amplitudes than the rest. 
These long-wave-length modes have velocity-space struc- 
ture close to Maxwellian, whose Hankel transform yields 
fJo(pv±)e-^ / ^v x dv x /v^ = e- p2, i /4 /2. The gradual 
steepening at the high-p region may be related to the gap 
in Fig.|6](Notice that Fig.|8]corresponds to taking the hori- 
zontal sum of Fig.|6]l. The wiggles in the high-p end of the 
runs C2 and D2 come from the slight lack of resolution, 
which goes away by increasing the velocity resolution as 
the high velocity-resolution run C3 shows. This increase 
of the velocity-space resolution does not affect the wave- 
number spectra, while failure to resolve the scaling regime 
does (see line CI in Fig.0. 

3.6 Entropy transfer 

In the scaling theory (Sec. 12. 2\ we assumed local-scale 
interaction by following Kolmogorov's argument. The ap- 
plicability of such an assumption may be directly investi- 
gated in the numerical simulation. 

In order to make this diagnostic we have introduced a 
Fourier filtered function (36l defined by 

ik-R 



gK(R,v x ):=Y,Sk(y x )e 11 



(20) 



ke'K 



where K = {k : Kp- 1 /2 < \k\p <Kp+l /2}. The filtering 
is orthogonal, and gK denotes the component of the distri- 
bution function with scale K l . Then the entropy amount 
contained in this scale obeys the evolution equation 

— f dRdv=Y Tf(K, Q) -collisions, (21) 
dtj 2F *g 

where Tf(K, Q) denotes the nonlinear transfer of entropy 
from scale Q l to K~ l 



Tf(K,Q) 



■-ff 



v±- 



gKV ExB ■ Vgg 



dRdv x , (22) 



which is by definition antisymmetric with respect to the 
exchange of two arguments. 

The numerical result obtained from the run D2 (see 
Table |2]i is shown in Fig. [9] which is normalized by 38 at 
each time and averaged over 10 < f/T; n j t < 15. Figure 
[9] shows a remarkable locality of the interaction, which is 
clearly seen even at each time. The color denotes the di- 
rection of the entropy transfer, which indicates that the en- 
tropy is transferred from large scales to small scales. 

The scaling theory does not require locality in velocity 
space scales, however the velocity-space transfer function 
may still give some useful information. It may be moni- 
tored by Hankel filtered function defined by 



gp(v x ) 
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Fig. 9 Time- averaged normalized entropy transfer function 
Tf(K, Q) [see also l[22]l] for the run D2 in Tabled 
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Fig. 10 Time-averaged normalized entropy transfer function 
T H (P,S) [see also CT l for the run C2 in Tabled 



where P = {p : Pv th - 1/2 < pv th < Pv lh + 1/2). Then, 
similarly to the Fourier filtering, entropy transfer function 
in Hankel space is described by 



T H (P,S):-- 



gpVExB ■ V| s v ± dv x dR, (24) 



gk(p)Mpv x )p dp, 



(23) 



where g(v x ) := g(v x )/ y/Fo is introduced to make entropy 
a bilinear form, and Th(P,S) = -Tn(S,P). 

Figure[10]shows the normalized, time-averaged trans- 
fer function in the velocity space. For numerical purposes 
the integral in d23l is approximated by a single representa- 
tive Hankel mode. Nevertheless the transfer function is lo- 
calized along the diagonal very well, which, together with 
Fig.|9l may suggest the entropy transfer along the diagonal 
in (k x , p)-p\ane (see Fig. [6]). 

4. Summary 

We presented electrostatic, decaying turbulence sim- 
ulations for weakly collisional, magnetized plasmas using 
the gyrokinetic model in 4D phase space (two position- 
space and two velocity-space dimensions; the extension to 
three spatial dimensions is left for future work). Landau 
damping was removed from the system by ignoring varia- 
tion along the background magnetic field. 

Nonlinear interactions introduce an amplitude- 



dependent perpendicular phase mixing of the gyrophase- 
independent part of the perturbed distribution function 
and create structure in v± which is finer for higher k± (see 
Fig. [6]). We found that the wave-number (Fourier) spectra 
of the perturbed distribution function and the resulting 
electrostatic fluctuations at sub-Larmor scales agreed well 
with theoretical predictions based on the interpretation 
of the nonlinear phase mixing as a cascade of entropy 
in phase space (see Figs. [5] and [7]) lfT2l l25l l27l . The 
velocity-space (Hankel) spectra show a rough consistency 
with the theoretical scaling, although the agreement is not 
as good as that of the wave-number spectra (see Figs. [8] 
and [(J. 

We introduced a dimensionless number (analogous to 
Reynolds number) that characterized the scale separation 
between the thermal Larmor scale and the collisional cut- 
off in phase space [see (TTTIl l. and showed that this number 
correctly predicted the resolution requirements for nonlin- 
ear gyrokinetic simulations (see Table|2]i. We also showed 
the trend that the entropy generation rate (or irreversibility) 
is independent of the collision frequency in the asymptotic 
limit of the weak collisionality (see Fig. [TJ. Finally we 
presented diagnostics of nonlinear transfer functions and 
showed that local-scale cascade of entropy is supported 
very well in both Fourier and Hankel spaces. 

We note that there are, in general, entropy cascades 
for each plasma species. Equations for the gyrokinetic tur- 
bulence at and below the electron Larmor scale are mathe- 
matically similar to the model simulated here and identical 
arguments apply 1 25 27 1. Similar considerations are also 
possible for ion-scale electromagnetic turbulence ll25l and 
for minority ion species (with some differences to be dis- 
cussed elsewhere). 

The structure of the small scales in phase space that 
we have presented is likely to be a universal feature of 
magnetized plasma turbulence. Understanding it theoret- 
ically and diagnosing it numerically is akin to the inertial- 
range studies for Kolmogorov turbulence, extended to the 
kinetic phase space. One should expect rich and interest- 
ing physics to emerge and it is likely that, just like in the 
case of fluid turbulence, predicting large-scale dynamics 
will require effective models for the small-scale cascade. 
An immediate physical implication of the existence of the 
entropy cascade is a turbulent heating rate independent of 
collisionality in weakly collisional plasmas. 
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